Degradation statistics
baRulho:quantifying habitat-induced degradation of (animal) acoustic signals
Source code, data and annotation protocol found at https://github.com/maRce10/barulho_paper
Purposes
Estimate degradation metrics on re-recorded signals from playback experiment at Bosque de Tlalpan, Mexico City, 2019
Run statistical analyses
Load packages
Code
source("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R")
# install/ load packages
load_packages(packages = c("remotes", "kableExtra", "knitr", "formatR",
"rprojroot", "maRce10/warbleR", "ggplot2", "tidyr", "viridis",
"corrplot", "brms", "ggdist", "cowplot", "cmdstanr", "maRce10/brmsish",
"emmeans", "ggsignif"))
my.viridis <- function(...) viridis(alpha = 0.5, begin = 0.3, end = 0.7,
...)
source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")Code
degrad_df <- read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")
degrad_params <- c("blur.ratio", "spectrum.blur.ratio", "envelope.correlation",
"excess.attenuation", "signal.to.noise.ratio", "cross.correlation",
"tail.to.signal.ratio", "tail.to.noise.ratio", "spectrum.correlation")
comp.cases <- complete.cases(degrad_df[, names(degrad_df) %in% degrad_params])
pca <- prcomp(degrad_df[comp.cases, names(degrad_df) %in% degrad_params],
scale. = TRUE)
# add to data
degrad_df$pc1 <- NA
degrad_df$pc1[comp.cases] <- pca$x[, 1]
degrad_df$pc1.1m.rate <- degrad_df$distance
# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:4])
pca_rot_stck <- stack(pca_rot)
pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)
ggplot(pca_rot_stck, aes(x = variable, y = rotation, fill = Sign)) +
geom_col() + coord_flip() + scale_fill_viridis_d(alpha = 0.7,
begin = 0.3, end = 0.8) + facet_wrap(~ind) + theme_classic()
1 Correlation among metrics
Raw metrics:
Code
cormat <- cor(degrad_df[, degrad_params], use = "pairwise.complete.obs")
rownames(cormat) <- colnames(cormat) <- degrad_params
cols_corr <- colorRampPalette(c("white", "white", viridis(4, direction = -1)))(10)
cp <- corrplot.mixed(cormat, tl.cex = 0.7, upper.col = cols_corr,
lower.col = cols_corr, order = "hclust", lower = "number", upper = "ellipse",
tl.col = "black")
2 Data description
- 20 frequencies
- 3 locations
- 11520 test sounds
- 160 sound treatment combinations
- Sample sizes per location, transect and signal type
Code
agg <- aggregate(cbind(sound.id, treatment.replicates) ~ location +
habitat.structure + distance, degrad_df, function(x) length(unique(x)))
agg$replicates <- agg$sound.id/agg$treatment.replicates
df1 <- knitr::kable(agg, row.names = FALSE, escape = FALSE, format = "html")
kable_styling(df1, bootstrap_options = c("striped", "hover", "condensed",
"responsive"), full_width = FALSE, font_size = 15)| location | habitat.structure | distance | sound.id | treatment.replicates | replicates |
|---|---|---|---|---|---|
| 1 | closed | 10 | 480 | 160 | 3 |
| 2 | closed | 10 | 480 | 160 | 3 |
| 3 | closed | 10 | 480 | 160 | 3 |
| 1 | open | 10 | 480 | 160 | 3 |
| 2 | open | 10 | 480 | 160 | 3 |
| 3 | open | 10 | 480 | 160 | 3 |
| 1 | closed | 30 | 480 | 160 | 3 |
| 2 | closed | 30 | 480 | 160 | 3 |
| 3 | closed | 30 | 480 | 160 | 3 |
| 1 | open | 30 | 480 | 160 | 3 |
| 2 | open | 30 | 480 | 160 | 3 |
| 3 | open | 30 | 480 | 160 | 3 |
| 1 | closed | 65 | 480 | 160 | 3 |
| 2 | closed | 65 | 480 | 160 | 3 |
| 3 | closed | 65 | 480 | 160 | 3 |
| 1 | open | 65 | 480 | 160 | 3 |
| 2 | open | 65 | 480 | 160 | 3 |
| 3 | open | 65 | 480 | 160 | 3 |
| 1 | closed | 100 | 480 | 160 | 3 |
| 2 | closed | 100 | 480 | 160 | 3 |
| 3 | closed | 100 | 480 | 160 | 3 |
| 1 | open | 100 | 480 | 160 | 3 |
| 2 | open | 100 | 480 | 160 | 3 |
| 3 | open | 100 | 480 | 160 | 3 |
3 Statistical analysis
Code
iter <- 10000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"), prior(cauchy(0, 4),
class = "sd"))
# set frequency to mean-centered
degrad_df$frequency <- degrad_df$frequency - mean(degrad_df$frequency)
# set base level for factors
degrad_df$habitat.structure <- factor(degrad_df$habitat.structure,
levels = c("open", "closed"))
degrad_df$frequency.modulation <- factor(degrad_df$frequency.modulation,
levels = c("no_fm", "fm"))
degrad_df$amplitude.modulation <- factor(degrad_df$amplitude.modulation,
levels = c("no_am", "am"))
degrad_df$duration <- factor(degrad_df$duration, levels = c("short",
"long"))
degrad_df$location <- as.factor(degrad_df$location)
degrad_df$distance_f <- paste0(degrad_df$distance, "m")
degrad_df$distance_f <- factor(degrad_df$distance_f, levels = c("10m",
"30m", "65m", "100m"), ordered = TRUE)
set.seed(123)
cmdstanr::set_cmdstan_path("~/Documentos/cmdstan/")
# to run within-chain parallelization
mod_pc1 <- brm(formula = pc1 ~ frequency * habitat.structure + frequency.modulation *
habitat.structure + amplitude.modulation * habitat.structure +
duration * habitat.structure + mo(distance_f) + (1 | location) +
(1 | treatment.replicates), data = degrad_df, prior = priors,
iter = iter, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
max_treedepth = 15), file = "./data/processed/regression_model_pc1.RDS",
file_refit = "always", backend = "cmdstanr", threads = threading(10))
mod_blurratio <- brm(formula = blur.ratio ~ frequency * habitat.structure +
frequency.modulation * habitat.structure + amplitude.modulation *
habitat.structure + duration * habitat.structure + mo(distance_f) +
(1 | location) + (1 | treatment.replicates), data = degrad_df,
prior = priors, iter = iter, chains = chains, cores = chains,
control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/regression_model_blur_ratio.RDS",
file_refit = "always", backend = "cmdstanr", threads = threading(10))3.1 PC1 degradation
Code
effects <- c("habitat structure", "frequency modulation", "amplitude modulation",
"frequency", "duration", "frequency:habitat structure", "habitat structure:duration",
"habitat structure:amplitude modulation", "habitat structure:frequency modulation")
mod <- readRDS("./data/processed/regression_model_pc1.RDS")
ggs <- extended_summary(mod, n.posterior = 1000, fill = viridis(10)[7],
trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
"frequency.modulationfm", "durationlong"), gsub.replacement = c("",
"habitat structure", "amplitude modulation", "frequency modulation",
"duration"), effects = effects, print.name = FALSE, trace = TRUE,
return = TRUE)3.1.1 Conditional plots
Code
bs <- 10
ts <- 4
ggf <- plot(conditional_effects(mod, "frequency:habitat.structure"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Open",
"Closed")) + theme_classic(base_size = bs) + labs(x = "Frequency",
y = "Degradation (PC1)", color = "Habitat") + scale_fill_discrete(guide = "none") +
theme(legend.position = c(0.2, 0.95)) + scale_x_continuous(breaks = seq(-5,
5, 2.5), labels = round(seq(-5, 5, 2.5) + mean(read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")$frequency),
1)) + geom_text(x = 0, y = -0.5, label = "*", size = ts)
ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Tonal",
"FM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "",
color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(c(-5, 1.5))
diff <- 0.1
# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(0), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggfm <- ggfm + geom_signif(y_position = c(1), xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
color = "gray40")
ggam <- plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Flat",
"AM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "Degradation (PC1)",
color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(c(-5, 1.5))
# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(0), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggam <- ggam + geom_signif(y_position = c(1), xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4,
color = "gray40")
ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Short",
"Long")) + theme_classic(base_size = bs) + labs(x = "Habitat",
y = "", color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(c(-5, 1.5))
# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(0), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("ns"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggdur <- ggdur + geom_signif(y_position = c(1), xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4,
color = "gray40")
plot_grid(ggf, ggfm, ggam, ggdur)
3.2 Blur ratio
Code
mod <- readRDS("./data/processed/regression_model_blur_ratio.RDS")
extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
"frequency.modulationfm", "durationlong"), gsub.replacement = c("",
"habitat structure", "amplitude modulation", "frequency modulation",
"duration"), effects = effects, print.name = FALSE)| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 6) Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) | blur.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) | 10000 | 4 | 1 | 5000 | 187 | 0 | 1772.987 | 2746.438 | 497531395 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| habitat structure | 0.029 | 0.026 | 0.033 | 1.001 | 14071.545 | 14362.080 |
| frequency modulation | 0.065 | 0.058 | 0.072 | 1.001 | 2201.419 | 3919.558 |
| amplitude modulation | 0.024 | 0.017 | 0.031 | 1.002 | 2032.491 | 4248.172 |
| frequency | 0.001 | 0.000 | 0.002 | 1.002 | 2355.250 | 4752.024 |
| duration | 0.001 | -0.006 | 0.008 | 1.002 | 1772.987 | 3901.099 |
| frequency:habitat structure | 0.003 | 0.002 | 0.003 | 1.001 | 7870.588 | 2746.438 |
| habitat structure:duration | 0.005 | 0.001 | 0.009 | 1 | 21697.716 | 14527.390 |
| habitat structure:amplitude modulation | 0.002 | -0.001 | 0.006 | 1.001 | 19309.028 | 14793.951 |
| habitat structure:frequency modulation | -0.020 | -0.024 | -0.017 | 1 | 22550.716 | 13648.971 |
3.2.1 Conditional plots
Code
bs <- 10
ts <- 4
ggf <- plot(conditional_effects(mod, "frequency:habitat.structure"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Open",
"Closed")) + theme_classic(base_size = bs) + labs(x = "Frequency",
y = "Degradation (PC1)", color = "Habitat") + scale_fill_discrete(guide = "none") +
theme(legend.position = c(0.2, 0.95)) + scale_x_continuous(breaks = seq(-5,
5, 2.5), labels = round(seq(-5, 5, 2.5) + mean(read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")$frequency),
1)) + geom_text(x = 0, y = 0.16, label = "ns", size = ts)
ylim <- c(-0.1, 0.5)
ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Tonal",
"FM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "",
color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
diff <- 0.1
# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(0.25), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggfm <- ggfm + geom_signif(y_position = c(0.3), xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
color = "gray40")
ggam <- plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Flat",
"AM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "Degradation (PC1)",
color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(0.2), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggam <- ggam + geom_signif(y_position = c(0.3), xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4,
color = "gray40")
ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Short",
"Long")) + theme_classic(base_size = bs) + labs(x = "Habitat",
y = "", color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(0.2), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("ns"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggdur <- ggdur + geom_signif(y_position = c(0.3), xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4,
color = "gray40")
plot_grid(ggf, ggfm, ggam, ggdur)
3.3 Excess attenuation
Code
mod <- readRDS("./data/processed/regression_model_excess_attenuation.RDS")
extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
"frequency.modulationfm", "durationlong"), gsub.replacement = c("",
"habitat structure", "amplitude modulation", "frequency modulation",
"duration"), effects = effects, print.name = FALSE)| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 6) Intercept-student_t(3, 83.8, 34) sd-cauchy(0, 4) sigma-student_t(3, 0, 34) simo-dirichlet(1) | excess.attenuation ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) | 10000 | 4 | 1 | 5000 | 9 | 0 | 2625.42 | 5199.726 | 865228054 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| habitat structure | 9.891 | 9.052 | 10.730 | 1 | 15473.413 | 14992.629 |
| frequency modulation | -3.158 | -4.709 | -1.594 | 1.002 | 2625.420 | 5199.726 |
| amplitude modulation | -2.724 | -4.204 | -1.227 | 1.001 | 3167.664 | 5593.286 |
| frequency | 3.164 | 2.902 | 3.431 | 1.002 | 2763.472 | 5433.162 |
| duration | -0.849 | -2.367 | 0.687 | 1.001 | 2700.931 | 5700.373 |
| frequency:habitat structure | -0.639 | -0.785 | -0.493 | 1 | 25778.804 | 13335.766 |
| habitat structure:duration | -0.399 | -1.241 | 0.427 | 1 | 20933.514 | 14843.771 |
| habitat structure:amplitude modulation | -1.576 | -2.414 | -0.739 | 1 | 21508.820 | 14345.505 |
| habitat structure:frequency modulation | -0.128 | -0.967 | 0.711 | 1 | 22068.512 | 14167.521 |
3.3.1 Conditional plots
Code
bs <- 10
ts <- 4
ggf <- plot(conditional_effects(mod, "frequency:habitat.structure"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Open",
"Closed")) + theme_classic(base_size = bs) + labs(x = "Frequency",
y = "Degradation (PC1)", color = "Habitat") + scale_fill_discrete(guide = "none") +
theme(legend.position = c(0.2, 0.95)) + scale_x_continuous(breaks = seq(-5,
5, 2.5), labels = round(seq(-5, 5, 2.5) + mean(read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")$frequency),
1)) + geom_text(x = 0, y = 70, label = "*", size = ts)
ylim <- c(20, 80)
ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Tonal",
"FM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "",
color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
diff <- 0.1
# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(65), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggfm <- ggfm + geom_signif(y_position = c(72), xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4,
color = "gray40")
ggam <- plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Flat",
"AM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "Degradation (PC1)",
color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(65), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggam <- ggam + geom_signif(y_position = c(72), xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
color = "gray40")
ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Short",
"Long")) + theme_classic(base_size = bs) + labs(x = "Habitat",
y = "", color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(65), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("ns"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggdur <- ggdur + geom_signif(y_position = c(72), xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4,
color = "gray40")
plot_grid(ggf, ggfm, ggam, ggdur)
3.4 Tail-to-signal ratio
Code
mod <- readRDS("./data/processed/regression_model_tail_to_signal_ratio.RDS")
extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE,
gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam",
"frequency.modulationfm", "durationlong"), gsub.replacement = c("",
"habitat structure", "amplitude modulation", "frequency modulation",
"duration"), effects = effects, print.name = FALSE)| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 6) Intercept-student_t(3, -6.3, 8.7) sd-cauchy(0, 4) sigma-student_t(3, 0, 8.7) simo-dirichlet(1) | tail.to.signal.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) | 10000 | 4 | 1 | 5000 | 16 | 0 | 2380.891 | 4713.849 | 1732769608 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| habitat structure | 3.253 | 2.897 | 3.605 | 1 | 13266.787 | 14830.753 |
| frequency modulation | -0.818 | -1.436 | -0.211 | 1.001 | 2503.012 | 5416.822 |
| amplitude modulation | 0.822 | 0.213 | 1.450 | 1.001 | 2720.556 | 5229.965 |
| frequency | 0.303 | 0.196 | 0.409 | 1.002 | 2749.731 | 4929.226 |
| duration | -0.426 | -1.029 | 0.182 | 1.001 | 2380.891 | 4713.849 |
| frequency:habitat structure | 0.020 | -0.041 | 0.082 | 1 | 25850.405 | 15259.070 |
| habitat structure:duration | -0.053 | -0.412 | 0.301 | 1 | 18539.674 | 15164.141 |
| habitat structure:amplitude modulation | -0.385 | -0.741 | -0.030 | 1 | 17231.026 | 14672.205 |
| habitat structure:frequency modulation | 0.540 | 0.183 | 0.898 | 1 | 20156.979 | 15407.009 |
3.4.1 Conditional plots
Code
bs <- 10
ts <- 4
ggf <- plot(conditional_effects(mod, "frequency:habitat.structure"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Open",
"Closed")) + theme_classic(base_size = bs) + labs(x = "Frequency",
y = "Degradation (PC1)", color = "Habitat") + scale_fill_discrete(guide = "none") +
theme(legend.position = c(0.2, 0.95)) + scale_x_continuous(breaks = seq(-5,
5, 2.5), labels = round(seq(-5, 5, 2.5) + mean(read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")$frequency),
1)) + geom_text(x = 0, y = 70, label = "*", size = ts)
ylim <- c(-27, -5)
ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Tonal",
"FM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "",
color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
diff <- 0.1
# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(-12), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggfm <- ggfm + geom_signif(y_position = c(-9), xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
color = "gray40")
ggam <- plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Flat",
"AM")) + theme_classic(base_size = bs) + labs(x = "Habitat", y = "Degradation (PC1)",
color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(-12), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("*"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggam <- ggam + geom_signif(y_position = c(-9), xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4,
color = "gray40")
ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"),
plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = c("Short",
"Long")) + theme_classic(base_size = bs) + labs(x = "Habitat",
y = "", color = "") + scale_fill_discrete(guide = "none") + theme(legend.position = c(0.93,
0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(-12), xmin = c(i -
diff), xmax = c(i + diff), annotation = c("ns"), tip_length = 0.02,
textsize = ts, size = 0.4, color = "gray40")
ggdur <- ggdur + geom_signif(y_position = c(-9), xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4,
color = "gray40")
plot_grid(ggf, ggfm, ggam, ggdur)
3.5 Combined results
Code
pc1 <- extended_summary(read.file = "./data/processed/regression_model_pc1.RDS",
n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
"habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
"durationlong"), gsub.replacement = c("", "habitat structure",
"amplitude modulation", "frequency modulation", "duration"),
effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)
gg_pc1 <- pc1$plot + labs(x = "Degradation (PC1)")
br <- extended_summary(read.file = "./data/processed/regression_model_blur_ratio.RDS",
n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
"habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
"durationlong"), gsub.replacement = c("", "habitat structure",
"amplitude modulation", "frequency modulation", "duration"),
effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)
gg_br <- br$plot + labs(x = "Blur ratio", y = "") + theme(axis.text.y = element_blank())
ea <- extended_summary(read.file = "./data/processed/regression_model_excess_attenuation.RDS",
n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
"habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
"durationlong"), gsub.replacement = c("", "habitat structure",
"amplitude modulation", "frequency modulation", "duration"),
effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)
gg_ea <- ea$plot + labs(x = "Excess attenuation", y = "") + theme(axis.text.y = element_blank())
tsr <- extended_summary(read.file = "./data/processed/regression_model_tail_to_signal_ratio.RDS",
n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis,
remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
"habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
"durationlong"), gsub.replacement = c("", "habitat structure",
"amplitude modulation", "frequency modulation", "duration"),
effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)
gg_tsr <- tsr$plot + labs(x = "Tail-to-signal ratio", y = "") + theme(axis.text.y = element_blank())
plot_grid(gg_pc1, gg_br, gg_ea, gg_tsr, nrow = 1, rel_widths = c(6,
2, 2, 2))
Code
estimates <- data.frame(pc1 = pc1$coef_table$Estimate, br = br$coef_table$Estimate,
ea = ea$coef_table$Estimate, tsr = tsr$coef_table$Estimate)
names(estimates) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
"Tail-to-signal ratio")
# make estimates relative to maximum estimate in data
rel_estimates <- data.frame(lapply(estimates, function(x) x/max(abs(x)) *
0.8))
names(rel_estimates) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
"Tail-to-signal ratio")
# corrplot(as.matrix(estimates), method = 'ellipse', )
signif <- data.frame(pc1 = pc1$coef_table$`l-95% CI` * pc1$coef_table$`u-95% CI` >
0, br = br$coef_table$`l-95% CI` * br$coef_table$`u-95% CI` >
0, ea = ea$coef_table$`l-95% CI` * ea$coef_table$`u-95% CI` >
0, tsr = tsr$coef_table$`l-95% CI` * tsr$coef_table$`u-95% CI` >
0)
names(signif) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation",
"Tail-to-signal ratio")
estimates <- cbind(rownames(pc1$coef_table), stack(estimates)[, 1],
stack(rel_estimates)[, 1], stack(signif))
names(estimates) <- c("predictor", "est", "relavite.est", "sig", "response")
estimates$relavite.est <- ifelse(estimates$sig, estimates$relavite.est,
0)
estimates$response <- factor(estimates$response, levels = rev(c("Degradation (PC1)",
"Blur ratio", "Excess attenuation", "Tail-to-signal ratio")))
estimates$predictor <- factor(estimates$predictor, levels = c("habitat structure",
"frequency modulation", "amplitude modulation", "frequency", "duration",
"frequency:habitat structure", "habitat structure:duration", "habitat structure:amplitude modulation",
"habitat structure:frequency modulation"))
ggplot(estimates, aes(x = predictor, y = response, fill = relavite.est)) +
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + geom_tile()
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + +
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + coord_equal()
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + +
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + #
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + use
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + blue
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + and
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + red
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + palette
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + scale_fill_gradient2(low
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + =
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + '#F36c68',
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + high
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + =
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + '#175493',
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + guide
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + =
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + 'none')
geom_tile() + coord_equal() + # use blue and red palette scale_fill_gradient2(low = '#F36c68', high = '#175493', guide = 'none') + +
scale_fill_gradient2(low = viridis(10)[3], high = viridis(10)[7],
guide = "none") + geom_text(aes(label = round(est, 2)), color = "black") +
labs(x = "", y = "") + theme_classic() + theme(axis.text.x = element_text(color = "black",
size = 11, angle = 30, vjust = 0.8, hjust = 0.8))
4 Takeaways
- Habitat structure seems to drive most of the observed degradation
- Frequency modulation and amplitude modulation were the signal structural features that more strongly affected transmission
- Signal features interact with habitat structure in diverse ways, sometimes increasing and sometimes decreasing degradation
- Degradation metrics are robust to variation in signal duration
5 Posterior predictive checks
Code
model_list <- c(pc1 = "./data/processed/regression_model_pc1.RDS",
blur_ratio = "./data/processed/regression_model_blur_ratio.RDS",
excess_attenuation = "./data/processed/regression_model_excess_attenuation.RDS",
tail_to_signal_ratio = "./data/processed/regression_model_tail_to_signal_ratio.RDS")
ndraws <- 20
for (i in seq_len(length(model_list))) {
cat("\n\n## ", names(model_list)[i], "\n\n")
fit <- readRDS(model_list[[i]])
ppc_dens <- pp_check(fit, ndraws = ndraws, type = "dens_overlay_grouped",
group = "habitat.structure") # shows dens_overlay plot by default
pp_mean <- pp_check(fit, type = "stat_grouped", stat = "mean",
group = "habitat.structure", ndraws = ndraws)
pp_scat <- pp_check(fit, type = "scatter_avg_grouped", group = "habitat.structure",
ndraws = ndraws)
pp_stat2 <- pp_check(fit, type = "stat_2d", ndraws = ndraws)
pp_loo <- pp_check(fit, type = "loo_pit_qq", ndraws = ndraws)
pp_error <- pp_check(fit, type = "error_scatter_avg_grouped",
ndraws = ndraws, group = "habitat.structure")
plot_l <- list(ppc_dens, pp_mean, pp_scat, pp_error, pp_stat2,
pp_loo)
plot_l <- lapply(plot_l, function(x) x + scale_color_viridis_d(begin = 0.1,
end = 0.8, alpha = 0.5) + scale_fill_viridis_d(begin = 0.1,
end = 0.8, alpha = 0.5) + theme_classic())
print(plot_grid(plotlist = plot_l, ncol = 2))
}5.1 pc1
5.2 blur_ratio
5.3 excess_attenuation
5.4 tail_to_signal_ratio
Session information
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
locale:
[1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
[3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
[5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
[7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ggsignif_0.6.2 emmeans_1.8.1-1 brmsish_1.0.0 cmdstanr_0.4.0
[5] cowplot_1.1.1 ggdist_3.2.0 brms_2.18.0 Rcpp_1.0.10
[9] corrplot_0.90 viridis_0.6.2 viridisLite_0.4.1 tidyr_1.1.3
[13] ggplot2_3.4.0 warbleR_1.1.28 NatureSounds_1.0.4 seewave_2.2.0
[17] tuneR_1.4.4 rprojroot_2.0.3 formatR_1.11 knitr_1.43
[21] kableExtra_1.3.4 remotes_2.4.2
loaded via a namespace (and not attached):
[1] backports_1.4.1 systemfonts_1.0.4 plyr_1.8.7
[4] igraph_1.4.3 splines_4.1.0 crosstalk_1.2.0
[7] TH.data_1.1-0 rstantools_2.2.0 inline_0.3.19
[10] digest_0.6.31 htmltools_0.5.5 fansi_1.0.4
[13] magrittr_2.0.3 checkmate_2.1.0 RcppParallel_5.1.5
[16] matrixStats_0.62.0 xts_0.12.2 sandwich_3.0-1
[19] svglite_2.1.0 prettyunits_1.1.1 colorspace_2.0-3
[22] signal_0.7-7 rvest_1.0.3 xfun_0.39
[25] dplyr_1.0.10 callr_3.7.3 crayon_1.5.2
[28] RCurl_1.98-1.12 jsonlite_1.8.5 lme4_1.1-27.1
[31] ape_5.6-2 survival_3.2-11 zoo_1.8-11
[34] glue_1.6.2 gtable_0.3.1 webshot_0.5.4
[37] distributional_0.3.1 pkgbuild_1.4.0 rstan_2.21.7
[40] abind_1.4-5 scales_1.2.1 mvtnorm_1.1-3
[43] DBI_1.1.3 miniUI_0.1.1.1 dtw_1.23-1
[46] xtable_1.8-4 proxy_0.4-27 stats4_4.1.0
[49] StanHeaders_2.21.0-7 DT_0.26 htmlwidgets_1.5.4
[52] httr_1.4.4 threejs_0.3.3 posterior_1.3.1
[55] ellipsis_0.3.2 pkgconfig_2.0.3 loo_2.4.1.9000
[58] farver_2.1.1 utf8_1.2.3 labeling_0.4.2
[61] tidyselect_1.2.0 rlang_1.1.1 reshape2_1.4.4
[64] later_1.3.0 munsell_0.5.0 tools_4.1.0
[67] cli_3.6.1 generics_0.1.3 ggridges_0.5.4
[70] evaluate_0.21 stringr_1.5.0 fastmap_1.1.1
[73] yaml_2.3.7 processx_3.8.1 purrr_1.0.0
[76] pbapply_1.7-0 nlme_3.1-152 mime_0.12
[79] projpred_2.0.2 xml2_1.3.3 brio_1.1.3
[82] compiler_4.1.0 bayesplot_1.9.0 shinythemes_1.2.0
[85] rstudioapi_0.14 gamm4_0.2-6 testthat_3.1.8
[88] tibble_3.2.1 stringi_1.7.12 ps_1.7.5
[91] Brobdingnag_1.2-9 lattice_0.20-44 Matrix_1.5-1
[94] nloptr_1.2.2.2 markdown_1.3 shinyjs_2.1.0
[97] fftw_1.0-7 tensorA_0.36.2 vctrs_0.6.2
[100] pillar_1.9.0 lifecycle_1.0.3 bridgesampling_1.1-2
[103] estimability_1.4.1 bitops_1.0-7 httpuv_1.6.6
[106] R6_2.5.1 promises_1.2.0.1 gridExtra_2.3
[109] codetools_0.2-18 boot_1.3-28 colourpicker_1.2.0
[112] MASS_7.3-54 gtools_3.9.3 assertthat_0.2.1
[115] rjson_0.2.21 withr_2.5.0 shinystan_2.6.0
[118] multcomp_1.4-17 mgcv_1.8-36 parallel_4.1.0
[121] grid_4.1.0 minqa_1.2.4 coda_0.19-4
[124] rmarkdown_2.20 shiny_1.7.3 base64enc_0.1-3
[127] dygraphs_1.1.1.6